home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / swagg_m.zip / MISC.SWG / 0125_The INFAMOUS Pentium Bug.pas < prev    next >
Pascal/Delphi Source File  |  1995-02-28  |  2KB  |  106 lines
  1. {
  2. I have been working with Cleve Moler (from MathWorks) and Tim Coe about a sw
  3. workaround for the Pentium FDIV bug.  Here is a BP version of the algorithm:
  4.  
  5. The FDIV_FIX unit defines a single function:
  6.  
  7.  
  8. function fdiv(x: extended; y: extended): extended;
  9.  
  10. which will use about 25% more cycles than a single fdiv call, but always return
  11. exact results for all double precision numbers, including the special cases of
  12. Nan, Inf and denormal.
  13. It also handles all extended precision numbers, giving a maximum error of a
  14. single bit in the last position (1ulp).  This error can only be introduced if
  15. the FDIV operations fails.
  16. During the unit startup code, I test the FDIV instruction, and if if fails,
  17. initialize a 16-byte table of critical nibble values.
  18.  
  19. If FDIV passes, the table is left empty, and no fixups will be done on
  20. divisions.
  21. =============== FDIV_FIX.PAS ===========================
  22. }
  23.  
  24.  {$N+,E-,G+}
  25.  Unit FDIV_FIX;
  26.  
  27.  Interface
  28.  
  29.  function fdiv(x, y : Extended): Extended;
  30.  
  31.  const
  32.    fdiv_risc_table : array [0..15] of byte = (0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);
  33.    fdiv_scale : single = 0.9375;
  34.    one_shl_63_l : Longint = $5f000000;
  35.  var
  36.    one_shl_63 : single absolute one_shl_63_l;
  37.  
  38.  Implementation
  39.  
  40.  function fdiv(x, y : Extended): Extended;
  41.  Assembler;
  42.  asm
  43.    fld [x]
  44.  @restart:
  45.    mov bx,word ptr [y+6]
  46.    fld [y]
  47.    add bx,bx
  48.     jnc @denormal
  49.  {$IFOPT G+}
  50.    shr bx,4
  51.  {$ELSE}
  52.    mov cl,4
  53.    shr bx,cl
  54.  {$ENDIF}
  55.    cmp bl,255
  56.     jne @ok
  57.  {$IFOPT G+}
  58.    shr bx,8
  59.  {$ELSE}
  60.    mov bl,bh
  61.    and bx,255
  62.  {$ENDIF}
  63.    cmp byte ptr fdiv_risc_table[bx],bh
  64.     jz @ok
  65.    fld [fdiv_scale]
  66.    fmul st(2),st
  67.    fmulp st(1),st
  68.     jmp @ok
  69.  @denormal:
  70.    or bx,word ptr [y]
  71.    or bx,word ptr [y+2]
  72.    or bx,word ptr [y+4]
  73.     jz @zero
  74.    fld [one_shl_63]
  75.    fmul st(2),st
  76.    fmulp st(1),st
  77.    fstp [y]
  78.     jmp @restart
  79.  @zero:
  80.  @ok:
  81.    fdivp st(1),st
  82.  end;
  83.  
  84.  const
  85.    a_bug : single = 4195835.0;
  86.    b_bug : single = 3145727.0;
  87.  
  88.  Procedure fdiv_init;
  89.  var
  90.    r : double;
  91.    i : Integer;
  92.  begin
  93.    r := a_bug / b_bug;
  94.    if a_bug - r * b_bug > 1.0 then begin
  95.      i := 1;
  96.      repeat
  97.        fdiv_risc_table[i] := i;
  98.        Inc(i,3);
  99.      until i >= 16;
  100.    end;
  101.  end;
  102.  
  103.  begin
  104.    fdiv_init;
  105.  end.
  106.